#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>
#include <time.h>


#define EPSILON 1.0e-5
#define REAL_T float
#define C 0.2

/* Pour obtenir le temps "horloge" : */
double my_gettimeofday(){
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return tv.tv_sec + (tv.tv_usec * 1.0e-6L);
}



int main(int argc, char **argv){
  int n = 0;
  int w = 0;
  int h = 0;
  int i = 0, j = 0;
  int step = 0;

  int grid_nb = 0, inv_grid_nb=0;
  REAL_T **grid[2] = {NULL, NULL};

  REAL_T max_error = 0.0;

  double t_debut = 0.0, t_fin = 0.0;
  double t_total_debut = 0.0, t_total_fin = 0.0;
  

  /* Valeurs par defaut : */
  w = 12;
  h = 300;

  /* Parametres : */
  if( argc > 1) h    = atoi(argv[1]);
  if( argc > 2) w    = atoi(argv[2]);

  /* Demarrage horloge : */
  t_total_debut = my_gettimeofday();

  /* Affichage informations diverses : */
  printf("Dims : (%d, %d) \t Nb points : %i \t Taille memoire : %f Mo X2\n",
	 w, h, w*h, (float) w*h*sizeof(REAL_T) / (1024*1024));

  /* Allocation memoire : */
  for (grid_nb=0; grid_nb<=1; grid_nb++){
    grid[grid_nb] = (REAL_T **) malloc(h*sizeof(REAL_T *));
    if (grid[grid_nb] == NULL){
      fprintf(stderr, "Echec allocation memoire pour 'grid[%d]'.\n", grid_nb);
    }
    for (i=0; i<h; i++){
      grid[grid_nb][i] = (REAL_T *) malloc(w*sizeof(REAL_T));
      if (grid[grid_nb][i] == NULL){
	fprintf(stderr, "Echec allocation memoire pour 'grid[%d][%d]'.\n", grid_nb, i);
      }
      
      /* set to 0: */
      memset(grid[grid_nb][i], 0, w*sizeof(REAL_T));
    } /* for i */
  } /* for grid_nb */

  /* Initialisation par defaut de la grille 0 : */
  grid_nb=0;
  for (i=0; i<h; i++){
    for (j=0; j<w; j++){
      grid[grid_nb][i][j] = (REAL_T) 25.0;
    }
  }

  /* Autres valeurs : 1 ligne de points a 100 toutes les 100 lignes */ 
  step = 100;
  for (i=step/2; i<h; i+=step){
    for (j=0; j<w; j++){
    	if(j==0)
	    	printf("grid[grid_nb][%d][%d] = 100.0;\n", i, j);
      grid[grid_nb][i][j] = 100.0;
    } /* for j */
  } /* for i */
  

  
   /* Test Exhaustif - Enregistrement des résultats */
	/*FILE * fi = fopen("results_seq_init", "w");
	 	 
	 for (i=0; i<h; i++){
	 	for(j=0;j<w;j++) {
	    	fprintf(fi, "%3.2f ", grid[grid_nb][i][j]);
   		}
   	
   		fprintf(fi, "\n");
    }
    
    fclose(fi);*/
  
  /***** Calcul : *****/
  do {
    inv_grid_nb = 1-grid_nb;

    /* Debut horloge pour cette iteration : */
    t_debut = my_gettimeofday();

    /* mise a zero de grid[inv_grid_nb]: */
    for (i=0; i<h; i++){
      memset(grid[inv_grid_nb][i], 0, w*sizeof(REAL_T));    
    } /* for i */

    for (i=0; i<h; i++){
      for (j=0; j<w; j++){
	
	grid[inv_grid_nb][i][j] = (1.0 - 4.0 * C) * grid[grid_nb][i][j]
	  + C * (j>0   ? grid[grid_nb][i][j-1] : grid[grid_nb][i][j])
	  + C * (j<w-1 ? grid[grid_nb][i][j+1] : grid[grid_nb][i][j])
	  + C * (i>0   ? grid[grid_nb][i-1][j] : grid[grid_nb][i][j]) 
	  + C * (i<h-1 ? grid[grid_nb][i+1][j] : grid[grid_nb][i][j]); 
	
      } /* for j */
    } /* for i */

    /* calcul de l'erreur */
    max_error = 0.0;
    for (i=0; i<h; i++){
      for (j=0; j<w; j++){
	REAL_T error = fabs(grid[grid_nb][i][j] - grid[inv_grid_nb][i][j]);
	if (error > max_error){ 
	  max_error = error; 
	}
      } /* for j */
    } /* for i */

    /* Arret horloge pour cette iteration : */
    t_fin = my_gettimeofday();

    /* Affichage pour cette iteration : */
    printf("Iteration %d : \t delta = %.3e \t temps = %.1f s\n", 
	   ++n, max_error, t_fin - t_debut);
	
    /* update 'grid_nb': */
    grid_nb = inv_grid_nb;
  } while (max_error > EPSILON && n < 40);
  /***** Fin du calcul *****/

	/* Test Exhaustif - Enregistrement des résultats */
	FILE * f = fopen("results_seq", "w");
	 	 
	 for (i=0; i<h; i++){
	 
	    	fprintf(f, "%3.5f ", grid[inv_grid_nb][i][1]);
   		
   		fprintf(f, "\n");
    }
    
    fclose(f);

  /* Liberation memoire : */
  for (grid_nb=0; grid_nb<=1; grid_nb++){
    for (i=0; i<h; i++){
      free(grid[grid_nb][i]);
    }
    free(grid[grid_nb]);
  }

  /* Arret horloge : */
  t_total_fin = my_gettimeofday();  

  /* Affichage des resultats : */
  printf("\nNombre d'iterations : %d\n", n); 
  printf("Erreur (delta) = %.3e\n", max_error); 
  printf("Temps total : %.1f s\n", t_total_fin - t_total_debut); 

  return EXIT_SUCCESS;
}


